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Abstract 

We calculate high-temperature graph expansions for the Ising spin glass model with 4 symmetric 
random distribution functions for its nearest neighbor interaction constants J^. Series for the 
Edwards- Anderson susceptibility Xea obtained to order 13 in the expansion variable { J /{k^T))'^ 
for the general d-dimensional hyper-cubic lattice, where the parameter J determines the width of 
the distributions. We explain in detail how the expansions are calculated. The analysis, using 
the Dlog-Pade approximation and the techniques known as Ml and M2, leads to estimates for the 
critical threshold {J/{k'QTc)Y and for the critical exponent 7 in dimensions 4, 5, 7 and 8 for all the 
distribution functions. In each dimension the values for 7 agree, within their uncertainty margins, 
with a common value for the different distributions, thus confirming universality. 

PACS: 05.70.Jk Critical point phenomena - 75.10.Nr Spin-glass and other random models 

1 Introduction 

In 1975 Edwards and Anderson (EA) [1] introduced a model for the theoretical study of spin glasses 
(SG) [2,3], which has started modern spin glass theory and has been of continued interest until today. 
Here we discuss the classical Ising case: The magnetic moments are represented by 'spin' variables 



1 



{si, i = 1,2, . . . , N}, which can assume the values Sj = ±1 and are located on the sites {i} of the d- 
dimensional hyper-cubic lattice. During our calculations we use a finite number of lattice sites N but 
eventually we are interested in the thermodynamic limit, — > oo. The spins' interaction is governed 
by the Hamiltonian 

AT 

'^{Jij}{Ui}) = - XI ^i^j -hoJ2^i^ (1) 

(ij) i=l 

where J2{ij) denotes the sum over all pairs of nearest neighbor lattice sites {ij), which we also call 
the lattice bonds, and the spin interaction constants Jij are chosen at random from a symmetric 
probability distribution, which is the same for all bonds. The external magnetic field is needed to 
define thermodynamic quantities as derivatives with respect to it, but apart from that, we concentrate 
on the case /iq = 0. The Hamiltonian's index {Jij} indicates that we deal with quenched disorder, 
i.e. the thermodynamic average for any observable is performed for a fixed set of coupling constants 
{Jij}- The configurational average of measurable thermodynamic quantities, over the random variables, 
is performed subsequently. For self averaging quantities this leads to expressions of what could be 
measured in experiments. We denote the thermodynamic average of any observable 74({si}) by {Ajx, 
and the configurational average of any function X{{Jij}) by [X]^^. 

The EA model neglects the details of the microscopic interaction between the spins, but exhibits 
the two essential ingredients that lead to the interesting features of spin glasses: Quenched disorder and 
frustration. Since little has been proved exactly for short ranged spin glass models, we assume what 
today is generally accepted, based on analytical and numerical evidence: Above the system's lower 
critical dimension di, whose value is controversial but agreed to be between 2 and 3 [4,5], it undergoes a 
continuous transition at a non-zero critical temperature Tc to a low temperature spin glass phase. This 
phase is characterized by broken spin-flip symmetry, i.e. a non-zero Edwards-Anderson order parameter 

1 ^ 

«- = ^E[(-^)t]r- (2) 

i=l 

The upper critical dimension, above which mean field behavior becomes dominant, is believed to be 
d« = 6 [6,7]. 

As the temperature T approaches Tc from above, we expect the susceptibility associated with q^^j^ , 
the Edwards-Anderson susceptibility, 

1 ^ 
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to exhibit a power law divergence, Xea ~ (^c — T)~^, characterized by the critical exponent 7. In the 
present study we use series expansions to investigate this behavior. Both q^^ and Xea ^'^^ related to 

configurational averages of higher order logarithmic derivatives of the partition function — ^q. 

""0 ho=o 

with respect to the external magnetic field. Those relations become linear in the thermodynamic 
limit [8]. 

The renormalization group theory [9] in dimension d = 6 — e predicts the universality of 7 and 
of other exponents, related to it by scaling relations. The universality classes should be set by the 
dimensionalities of space and of the spin variables, and not by details of the distribution functions. 

Series expansion has been used in the past to study the spin glass transition [8, 10-13] and the 
results support the statements mentioned above. Our renewed interest in the problem awoke with a 
series of studies [14-20] that found, based on computer simulations, that the critical exponents vary 
with the probability distribution for the quenched disorder in the coupling constants Jij. This is in 
clear violation of universality and not sufficiently explained by theory. 

Undoubtedly, many of the enormous complications and features observed in the study of spin 
glasses arise from the disorder inherent in these systems. They gave the model the reputation of being 
one of the toughest subjects in computational physics. Simulations are here directly impacted by long 
relaxation times, memory effects, hysteresis, the rugged energy landscape with many meta-stable states 
and the huge parameter space over which to average. 

The technique of series expansion comes with two immediate advantages: The averaging over the 
randomness can be done exactly, and the series can, given the availability of graph data, be obtained 
in general dimension. The subsequent analysis is still done in each dimension separately, but results 
generally get more reliable with increasing dimension, while simulations become increasingly expensive 
in their computational demands. The previous series expansion studies of the Ising spin glass used 
only the bimodal random distribution of Jij = itJ, limiting their use in the comparison with the 
claims of violated universality. In the present study we extend the research by addressing several other 
symmetric distribution functions, each with a variable width determined by the parameter J. We 
use the same distributions as Bernardi and Campbell in [15], except for the exponential distribution, 
which is excluded for reasons given in Sec. IHl After introducing additional notations and the random 
distribution functions in Sec. [21 we give a detailed explanation of the series generation in Sees. and 
^ which should allow the interested reader to follow each step. As an example, we actually show the 
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complete calculation of a fourth order series in Sec. [S| In Sec. IHl we present our general-dimension 
series in full, accompanied by some discussion of accuracy checks. Our series analysis and final results 
are described in Sec. [7| and we finish with our conclusions in Sec. |HJ 



2 Further Notations and Definitions 

With P = where fee denotes Boltzmann's constant and T the absolute temperature, the ensemble 
average of an observable A is calculated by 

= -V-^ = -fr(F«T- M 

where the partition function Z appears in the denominator. Here the trace (Tr) is a shorthand for 
summing over all possible values of the spins' {si} configuration 

TrX= Tr X({.a)= ' ' ' H ' (5) 

The free energy per site F is obtained from Z by 

F = —Fn = —\lnZ\ii. (6) 

Since the interaction constants Jij appear only in products with /3, it is convenient to use Kij = (3Jij 
as the argument of the distribution functions introduced below. If is some measure of [Jj^]r, then 
we also use k = (3 J as expansion variable, at least temporarily. Since only even powers of k remain, 
we eventually use x = k? as the expansion variable in our high temperature series. Likewise we use 
= {J/{kBTc))'^ to denote the critical threshold. 

In the general case of a continuous probability distribution -P(x), the configurational average is the 
nested integral 

— CO 

[X]r = /•••/ X{{x,,}) l[{P{x,,) dx,,) . (7) 

-oo fo'> 

For the bimodal random distribution the coupling constants Kij for nearest neighbor pairs randomly 
assume only values of either +k or — so the latter integral can be written as the nested sum 

[^]R = i E Xi{Ki,}), (8) 
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where a normalization factor of 1/2 stems from each Hij in the sum. In the d-dimensional hyper-cubic 
lattice with N sites the number of nearest neighbor pairs approaches Nd for large N when boundary 
effects become negligible. 

Near the critical temperature Tc, the quantity of our interest, Xea^ expected to have a singularity 
of the form 

Xea ~ ^(^c - xr^il + B(xc - + ■••)• (9) 

The aim of our analysis is to determine the critical exponent 7 and, to a lesser extent, the first correction 
exponent Ai. As for the free energy, we study this susceptibility per lattice site. 

2.1 The Random Distributions 

The different probability distributions, that we study, are listed below. We call them Bimodal, Gaus- 
sian, Uniform and Double-Triangular. 



{6[z - k] + 5[z + k]) 



^'bim(^) = 
-fgau {z) = 

Puni{z) = 

• The distributions are largely characterized by their moments 

/oo 
z^ P{z) dz. 
-00 




(10) 
(11) 
(12) 

(13) 



(14) 



• Since all distributions have the symmetry P[—z) = P{z), the moments for odd n vanish. In 
particular, the distributions have zero mean [z]r = z P{z) dz = 0. 

• For even n the moments are: 

Bimodal distribution: 
Gaussian distribution: (n — 1)!! 
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Figure 1: The distribution functions for k = 5. 
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' Double-Triangular 



Uniform distribution: 

Double- Triangular distribution: K"/(n/2-|-l) 

• Thus ah the distributions are normalized f^^P{z) dz = 1. 

• A distribution's second moment M2 = z'^ P{z) dz, equal to the variance, is commonly 
associated with its width. In all cases it is linear in k^, but with different pre-factors. Explicitly, 
M2 is equal to (Bimodal), (Gaussian), k^/S (Uniform), and (Double- Triangular) , 
respectively. With slightly redefined variables, M2 could be equal to in all cases, which, in 
retrospective, would have been nicer. 

• Figure n illustrates the distribution functions. The plots were calculated for the parameter k = 5. 



2.2 Tangential Moments 

In the calculation of the series we encounter the following integrals over the distributions 

nin = [/c"]r = j tanh"(z) P{z) dz. (15) 

We sometimes refer to m„ as the n-th tangential moment of the distribution, in order to distinguish it 
from the regular moment ()14() . For a series to power we need all moments up to m/vr (not 2N as 
we will see later). For the bimodal distribution the tangential moments are trivial: 

1 I odd n 

rUn = -(tanh"(K) + tanh"'(-K)) = < 

I tanh" K even n. 
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The simple form m2n = tanh^" k = w"' makes w = tanh^(K) an alternative (and convenient) expansion 
variable for this case, which has been used in the past [8, 10]. This, however, is not true for the other 
distributions. For them it may be possible to calculate the m„ analytically, as well. But the results may 
be complicated functions of the parameter k, not suitable for our power series expansion. Hence we are 
content to calculate the necessary moments m„ of all distributions as series in k. To avoid the tedious 
work, this can conveniently be done with software for symbolic computation, such as Mathematica. 
The obtained coefficients are later used in the graph expansion, both for the computerized calculation 
and for the example in this article. The symmetry of the distributions makes all moments for odd 
n vanish. Due to the power series of tanh" k, each moment m„ has only powers of and higher, 
somewhat important during cumulant subtraction. Note that in the framework of series expansions, 
this is an exact treatment of the randomness. We do not lose any additional information since a priori 
we are limited to the highest order of our final series. 

As an illustration, wc show the expansion of the first few moments for the bimodal and for the 
Gaussian distributions, to be used in the example below. 
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Bimodal Tangential Moments 



mo = i , 

17/^6 62 Ac^ 1382^10 _ 21844^12 

- ~3~^^5 315~^ 14175 467775 ^ 

929569 K^'^ _ 6404582 443861162 18888466084^20 

42567525 638512875 ^ 97692469875 9280784638125 
113927491862^22 58870668456604^24 8374643517010684 «;26 

+ 



126109485376875 147926426347074375 48076088562799171875 

689005380505609448 k^^ 129848163681107301953 
9086380738369043484375 ^ 3952575621190533915703125 ' 
4 4k6 6k8 848At^o 8507^12 

= "-^ + — -^4^ + ^417^- 

3868 K^'' 46471426 ^ 47060768 518299498^20 _ 

10395 212837625 383107725 ^ 7753370625 ~ 
92014385608^22 39319617599924^24 12160377940064^26 

2598619698675 ^ 2143861251406875 1304465840803125 
14121349128787207129 k^^ 20894145609081223868 k'^o 



3028793579456347828125 9086380738369043484375 ' 
_ 6 98 2266^12 1901^1'' 79214^16 136750052^^8 

me - K - K +^ 945 ^ 945 51975 ^ 127702575 

64742312^20 3282022 «;22 710423622556^24 82292419438259 «;26 

91216125 7309575 2598619698675 ^ 510443155096875 

68433004067940682 ^28 157107220075270779857 



739632131735371875 3028793579456347828125 



2.2.2 Gaussian Tangential Moments 



mo = 1 , (20) 

n A 17k^ 62 1382^10 21844^12 

1712 = K - 2 H \ h 

3 3 15 45 

929569^1'' 6404582^16 443861162 18888466084 k^o 

315 315 2835 14175 

1936767361654 58870668456604 8374643517010684 

155925 467775 ^ 6081075 

689005380505609448^28 129848163681107301953^3° 



42567525 638512875 

93577 fi;i2 



rrn = 3k''-20k^ + 126 - 848 + 
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14 46471426^16 800033056 ^i^ 9847690462^2° 

50284 ^1^ + + 

105 189 225 

92014385608 904351204798252 k^^ 12160377940064 ^26 

189 ^ 155925 165 

14121349128787207129 ^28 605930222680755492172 k^o 



14189175 42567525 

me = 15 210 + 2331^1° - 24926 + 271843 K^^- 

ifi 2324750884^18 20911766776^2° 

3089346^16 + + 

63 45 

^iTo.ooooo 22 16339743318788 k24 1892725647079957^26 

617o483o82 k 1 

189 1485 

68433004067940682 ^28 4556109382182852615853 k^o 



(22) 



+ ■ (23) 



3465 14189175 



3 Connected Graph Expansion and Cumulant Subtraction 

An extensive physical quantity X c an be expanded in terms of connected graphs only [21]. To order 
n in a suitable expansion variable, say x, all connected graphs with n or less edges are used, 

n 

^ = EE^(r)x^+o(ixri). (24) 

Here Ti, denotes the set of all connected graphs T with b edges, w{T) is the lattice constant for weak 
embeddings of T inside the lattice (see Sec. 13. 2j) and is the cumulant of the graph's contribution. 
For the quantity we calculate, the empty graph and the single vertex graph (with edges) together 
only contribute a constant summand of 1 with the chosen normalization. 
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Figure 2: The smallest graphs, that can be embedded into the hyper-cubic lattice. 
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The cumulant contribution of a connected graph T is obtained by subtracting off the cumulant 
contribution of all its connected subgraphs, 

X^ = Xr-^X^^. (25) 

Due to the subtractions, is the contribution to Xr, which depends on every one of the b edges in 
r and thus has only terms of order b and higher in the expansion variable. This property allows us to 
stop the expansion at a certain size of graphs, with a series which is correct to the obtained order, and 
is in contrast to the original Xr, which can contribute to any power. 

We often use the term bond instead of edge and likewise site instead of vertex since we deal with 
a physical model on a lattice, and will eventually embed the graph inside it. The physical model is 
also the reason we do not address digraphs or graphs with loops. Here the term loop, as commonly 
used in graph-theory, denotes an edge whose both ends are incident on the same vertex. This must be 
distinguished from a cycle (closed path) in the graph, which is important to us. 

Our model involves only nearest-neighbor interactions, visualized by occupied lattice bonds. Since 
a graph T on the lattice is completely isolated from the rest of the infinite lattice by unoccupied bonds, 
the thermodynamics of its spins is determined by the reduced A'^r-particle Hamiltonian for the graph 

nr{si G r} = - J^J^i^J- (26) 

3.1 Cumulant Subtraction 

Equation ()25() contains the sum over connected subgraphs. In the following we have written out the 
cumulants for the smallest graphs, with explicit numerical coefficients and graph indices, because we 
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will need them in the example later. For these small graphs, the expressions can easily be confirmed 
by visual inspection using Fig. [3 These cumulants are given by 





= ^To 












= Xr, 












= Xr^ 
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-2X^^ 
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- 3X^2 
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-4X^^ 




-^x^u 




= Xr, 


- 5^ro 


-^X^r. 




-^x^u 




= Xr^ 


- 5^ro 


-4X^^^ 




-1X^3 




= Xr^ 


-5^ro 


-4Xf^ 


-6^r. 


-4X^3 



(27) 
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3.2 The Lattice Constants for Weak Embeddings 

In (|24() . the cumulant contribution of each graph is multiplied by its lattice constant w{T). This 
constant is the number of distinct ways per lattice-site in which the graph can be weakly embedded in 
a particular lattice, and thus it ties our series to that specific lattice. In this article we address the 
d-dimensional hyper-cubic lattice, which, as the term suggests, is a generalization of the square lattice 
{d = 2) and the cubic lattice (d = 3). We use the tabulated functions wr{d) from [22]. To calculate 
our example later on, we need the lattice constants of the first graphs: 



w{To 
w{T2 
w{T4 

w{Tj 
w{T8 



1, 
d, 



4d3 
3 

d-Ad'^+Ad^ 



f -2^2+ 3 , 



(28) 



_d , ^ 
2 2 ' 

d + 4d2 _ I2d^ + 8d^, 
-2 d+ 10^2 - 16d^ + 8d^, 
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4 General Calculation and Simplifications for the Ising Spin Glass 

The Boltzmann factor e~^^ can be rewritten in a way that makes the calculation of the trace more 
convenient. With the notations A = Pho, r = tanhA, Kij = cosh Kij and kij = tanh we obtain: 

/ 

e = ex.p I P^JijSiSj + Php Sj 

N 

= (cosh(KjjSjSj) + sinh(KjjSjSj)) JJ (cosh(Asj) + sinh(Asi)) . (29) 

Now we exploit the fact that the only possible values of Si are ±1, together with the symmetry of cosh 
and anti-symmetry of sinh: 

N 

e-/3W = cosh^ A H Kij J] (1 + SiSjkij) J] + ^i'^) • (^0) 

{ij) {ij) «=1 

We address the case without external magnetic field, ho = 0, and call 

(31) 



= ^ Tr ( n (1 + kijSiSj) n (1 + SiT) ) 



the reduced partition function and sometimes use the notation ZR(r„) = Z^^^. For the case of zero 
external magnetic field /iq = we end up with the important (for our coming calculations) equation 



{A)T = 7 (32) 

^{s,} [U{ij) (l + kijSiSj) 



Often A will be a linear combination of SmSn- We will then see the terms 

2"^ Tr|3j (smSnWuj) (1 + hjSiSj) 



. y -l(y> „ „ , 

{SmSn)T = 7 ^ ■ (33) 

2-^Tr{,^} (n(.,) {l + h,s,Sj) 
4.1 Graph- Expansion for the Edwards- Anderson Susceptibihty 

For the graph expansion of the Edwards- Anderson susceptibility Q we apply the general formula ()24() . 
The calculation is done to order n in x = k^. With this expansion variable all the dependence on the 
coupling strength parameter J and the temperature T are absorbed in the argument of the power 
series. We obtain 

n 

^Xea = E E ^(r)Xf;. (34) 
12 



and also denote the associated (non-cumulant) observable on only one (sub)grapli as 

Xr=Yl [(««^j)t]r- (35) 

5 Explicit Calculation for the Smallest Graphs 

We now show explicitly the calculation of Xea smallest graphs in Fig. |21 using the equations 

from the previous section. Here we often use /c"- = following definition (jlSf) . 

L •'J R 

• For Fq, we have 

X,jFo) = l[(s2)2.]^=[(l)2]^ = l. (36) 
For any spin, the self-correlation trivially equals 1. 

• For Fi, 

XEA(ri) = [{siS2fT]^+ [(s2Si)|]r + {Slf + {S2f = 2 [(siS2>|]r + 2. (37) 

The denominator of (siS2)t is 

Zr, = Zjj(Fi) = 2-2 Tr (1 + fcississ) = 1, (38) 

{si,S2} 

and the numerator is 

(siS2)t Zr, = 2-2 Tr {siS2{l + /C12S1S2)) = 2-^ Tr {siS2 + ku) = ku- (39) 

{si,S2} {si,S2} 

Thus, 

X^A^i) = 2 [{siS2fT]^ + 2 = 2 [kl2\^ + 2 = 2 [tanh2(/3Ji2)]j^ + 2 = 2^2 + 2, (40) 
where we used the definition of the tangential moment 7712. 

• For F2, 

XeJF2) =2([(siS2)|^]r+ [(Wt]r+ [(^1S3)t]r) +3. (41) 
The denominator of each correlation is 

Zr, = 2-3 TV {{1 + kusis2)il + k23S2S3)) (42) 

{si,S2,S3} 

= 2-3 Tr (1 + /C12S1S2 + A:23S2S3 + ki2k23SiS3) = 1. (43) 

{si,S2,S3} 
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The numerator is 

{siS2)t Zr^ = 2"^ Tr (siS2(l + A;i2SiS2)(l + /C23S2S3)) = ku, (44) 

{si,S2,S3} 

and likewise 

{s2S3)t Zr^ = k23 , and (45) 
(siS3)t Zr^ = 2"^ TV (siS3(l + A;i2SiS2)(l + ^235253)) = fei2A;23- (46) 

{si,S2,S3} 

The configurational average gives 

[(siS2)t]r = [/c?2]r = "^2 , [(s2S3)t]r = m2 and (47) 

[{siss)l] J, = [kl,kl,] ^ = [kl,] ^ [kl,] ^ = ml, (48) 

where in the last step we used the fact that the random variables kij for different bonds are 
uncorrected . Finally, 

XEA(r2) =3 + 4m2 + 2mi. (49) 

A few general observations are in order: Prom the possible values of a spin variable, = ±1, we trivially 
have s| = 1. Further, the trace sums over the possible values of each spin, so any summand with an 
unpaired Si vanishes. When writing out the product within the traces (1 + kijSiSj) •••(! + kkis^si) one 
has to choose from each pair of (ij) (and for any resulting summand) either the constant 1 or the bond 
factor kij. To find the terms that will actually survive the trace, the SiSj site-factors accompanying 
any bond have to be combined with other sites, either from another bond in the graph or from terms 
originally present inside the trace. 

For any acyclic graph, which we have seen so far, we immediately see that the reduced partition 
function is always equal to 1; no bond can be combined with others to eliminate all unpaired sites. Also 
the numerator for pair correlations remains simple. Visually only a path of bonds joining the two sites 
in question yields factors of only paired spins, which survive the trace. Each such bond contributes a 
factor of m2. A constant multiplier results from the number of ways the pair of sites (and equivalent 
pairs) can be joined. For cyclic graphs the calculation becomes much more complicated. In passing 
please note that one of the inherent features of spin glasses is frustration. Only a cyclic graph is 
susceptible to this phenomenon and as such can possibly integrate true SG-properties into our series. 

14 



For Ta, 



XEA(r3) = 2 [(siS2)t + (s2S3)| + (s2S4)t + ('SiS3)t + (si54)t + ('53S4)t]r + 4 

= 4 + 6m2 + Gmn. 



(50) 



For 



XEA(r4 



2( [(siS2)|]r+ [(s2S3)|]r+ [(s3S4>r]R + 

[(siS3)|]r + [(s254>t]r + [(siS4>r]R ) + 4 

4 + 6m2 + 4m2 + 2m2. 



(51) 



• For Fs, 

X^jFs) =4 + 8 [(.^153)1] J, + 4 [(siS2)|]r . (52) 

Here we used the equivalence of pairs of sites, in terms of adjacencies, to reduce the number of 
terms to 3. We defer calculation of the result till later. 

• For Fe, 

X^aO^g) = 2 [4(siS4)^ + 3(siS2)| + 2(siS5)| + l{siS3)l]^+5 = 5+8m2+6mi+4mi+2m^. (53) 



For Ft, 



XeaO^y) = 2 [4(siS4)t + 4(siS2)t + 2{siS5)'^]^ + 5 = 5 + 8m2 + 8mi + 4ml (54) 



For Fr 



XEA(r8) = 2 [4(siS5)^ + 6(siS2)t]r + 5 = 5 + 8m2 + 12m2. (55) 



We now resume the calculation of Xea(-'^5) from which, due to the graph's cycle, is significantly 
more complicated than the contribution of the other graphs. The complexity enters through the non- 
trivial partition function 



Zr, = 2 



-4 



Tr ((1 + A;i3SiS3)(l + A;23S2S3)(1 + kuSiSi){l + ^245254)) = 1 + ki3k23kuk2i. (56) 

{si,S2,S3,S4} 

The numerators for the spin correlations are {siSs)t Zr^ = ^13 + fci4A;24^23 and {siS2)t = ki3k23 + 
^14^24- The complicated part is now performing the averaging over the randomness for 



[(siS3>r]R 



fcl3 + ki4k24k23 
1 + ^13^23^14^:24 
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(57) 



R 



and 

" - ^' - - I 

1 + ki3k23ki4k24 ) 



2" 

^13^23 + ^14^24 \ 



(58) 

R 

In fact, for most probability distributions we do not know to calculate this directly. For the bimodal 



distribution it is possible, and was in fact done also for much larger graphs with a computer [8]. For 
continuous distributions the calculation would at best become extremely tedious and most likely not 
feasible for large graphs. Our solution again utilizes a power expansion. By this we do not lose any 
more information since our final series are limited to a certain order in the expansion variable anyway. 

To make the expansion process more obvious we rewrite the equations with the symbols uij = kij/k 
such that kij = k Uij and 

1 + k^Ui3U23y'14U24 ' 

/ V _ k'^Ui3U23 + fc^ni4tt24 /g^x 
1 + k'^Ui3U23Ul4U24 ' 

and expand the squares in powers of k: 

(siS3)t = ul^k^ + 2ui3UuU23U24k^ + i—'2Ul3UUU23U24 + 'Wl4'"23^24)'^^ 
-4(^^?3"l4«23^24)^^ + {Suj^uj^ul^ul^ - 2«i3U?4«^3ui4)fc^° 

+6ul^ul4ul^ul4k^^ + {-Aul^ul^ul^uli + Sul^uf^ui^ui^jk^^ + 0(|A;|^^), (61) 

{siS2)t = k^{ulsul^ + 2U13UUU23U24 + UuuIa) 

+A;^(-2n?3ni4U23ti24 - 4:ul^ul4^ul^ul4^ - 2ni3uf4n23«i4) 

+k^\3uisul443ul4 + Gul^ul^ul^ul^ + Sul^uf^ulsu^^) + 0(|A;|^^). (62) 

With the fractions removed, we can again factorize for averages over independent variables, and thus 
use the previously defined moments of the random distributions: 

[(*i«3)t]r = m2 + m| - + 3m|m4 + 3m2m| + 0(|k|^^), (63) 
[(siS2)t]r = 2ml-4mi + 6mlml + 0{\K\^^). (64) 

Here we expanded to a higher power than actually necessary for the largest graph that we consider in 
this example. It shows that higher moments actually show up. We first encountered 1113 = and 
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For the bimodal distribution we can quickly use the simple form of the moments m2 = w and m4 = nP' , 
and obtain 

[{sis^)t] r = w + w'^ -Aw^ + Zw^ + 2>w^ + ... (66) 
[(siS2)t]r = 2u;2 - 4u;^ + 6t(;'5 + . . . , (67) 

and thus 

X^p^{T^^) = 4 + S-u; + Su-^ + - 48u'^ + 24w;^ + 24u;'^ + 24u;^ + 0{\w\^). (68) 
For the general distribution we remain with 

^^^^(rs) = 4 + 8(m2 + ml + - 6m2 + 3m^m4 + Smgrnl + 2,m2m\) + 0{\k\^^), (69) 

on which we elaborate further. 

5.1 Performing Cumulant Subtraction 

In this section we perform the cumulant subtraction for the quantity Xea' which we calculated in 
the previous section. The general equations were given in Sec. 13.11 and in the following we show the 
cumulant graph contributions first in terms of the tangential moments and also for the case of the 
bimodal distribution using w as the expansion variable, thus substituting m2 = w and m4 = w"^. 





= 1, 




= 2m2 = 2w, 




= 2ml = 




= 0, 




= 2ml — '^w^i 




= —487712 + 2Amlmi + 2Amlm\ + 2Am2m\ + 




= -ASw^ + 2Aw^ + 2Aw^ + 2Aw'^ + . . . , 




= 2ml = 2'"^^' 




= 0, 




= 0. 



(70) 



Indeed we see that no graph contributes to a power of w less than its number of bonds. This fact is 
used in the computerized calculation as an internal check. 
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In general we use k (or x = k^) as the expansion variable, for which we now use the expansions of 
the moments nin from section ^ITR We show the result for the Gaussian distribution: 

X^^ = 2k4_8k6 + 5^ + O(|k|10), 

= 0' 

X^^ = 2^6- 12^8 + O(|k|10), (71) 

= -48^8 + O(|k|10), 

X^^ = 2k8 + O(|k|10), 

^rs = 0. 
5.2 Using the Lattice Constants 

Using ()24() together with the lattice constants of the smallest graphs, as given in Sec. 13.21 we can now 
perform the final summation and obtain the series. For the Gaussian distribution, we find 



+ ( _58d+^^- 72^3 + 16^4 + (72) 



For the bimodal distribution one has 



= l + 2dK^ +i^- — + Ad^)K''+{— — + 8^3)^^ 

/1210d , 24d2 ,„ . , i«^4\ 8 



+ — — + ^-40d^ + 16dMK« + ... (73) 
\ d3 5 J 

or, expanded in w, 

= 1 + 2 du- + (-2 d + 4 d^) + (2 d - 8 + 8 d^) u;^ + (26 d - 16 - 24 + 16 d^) u-^ + . . . (74) 



6 The Full Series 



The series for the d-dimensional hyper-cubic lattice to order 13 need to take into account 20724 graphs 
of up to 13 edges, and are hence calculated using computers. We use the graph data files that were 
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originally prepared for [22] by Wan et al, and have since been used in many studies. Programs were 
written, that use these data files to compute the series as outlined in the previous sections. Details of 
the algorithms, including important efficiency considerations, are presented in [23]. 

TablesHto^show the resulting series in full. In Tab.Q for the bimodal distribution, the coefficients 
are given as exact fractions. For the other distributions they were in part calculated using the data 
type long double in C++ which limits their accuracy. For comparison we used two different processor 
architectures where this data type is represented in either 96 or 128 bits, and also compared part 
of the data with calculations done in double precision (64 bits). Small rounding errors are obvious 
in most numbers, but further investigation shows, that for coefficients large in absolute value, the 
numerical accuracies become important. Originally we had included in this work the exponential 
random distribution which is also addressed in [15]. This distribution decays slower than the others 
and the resulting coefficients become very large in absolute value, to a degree that intermediate numbers 
either can not be presented in long double variables or the rounding errors become so dominant that 
the highest order coefficients come out completely wrong. We have started to calculate the series using 
arbitrary-precision numerical libraries, but that work was not ready in time to be included here. We 
exclude the exponential distribution from the present work, and for the remaining series present the 
coefficients in as many digits as we expect to be correct from the comparisons mentioned above. 

From experience we know that small changes in the coefficients do not influence the results obtained 
from series analysis. Hence the numerical inaccuracies present in the power series should not influence 
our final results. For the bimodal distribution we supplement coefficients for orders x^'* and x^^, 
which were calculated using the non-free-end (NFE) technique and associated graph data. In this 
technique by Harris [24] the thermodynamic functions under study are renormalized in such a way 
that the contribution from a graph with at least one free end (i.e. a vertex with only one incident 
edge) vanishes. This renormalization is possible for the bimodal distribution [8] but was not obtained 
for the others. We use equations from [8] for the NFE-expansion in w, but do not describe the process 
here since the series can also be obtained directly by variable transformation from w to x, which indeed 
we use as a consistency check. 

Several checks are performed to assure the correctness of our series expansions: The first is a 
complete recalculation of the corresponding series in [8] for the bimodal distribution, which shows that 
the algorithm and its implementation are basically correct. 
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We mentioned earlier that after cumulant subtraction, a graph of b edges has only terms of order 
b and higher in the expansion variable. As an additional check we do the actual calculation of the 
vanishing terms, track the maximal deviation from zero, and confirm that this number is in the same 
range as the numerical rounding errors observed elsewhere. 

For a few sequences of coefHcients we find, by examination of the numerical values, what their exact 
value must be in general. If we denote by aij the coefficient multiplying x'^dP we observe: 

• For the bimodal distribution an = 2* and aj = —(5/6) 2* {i — 1). 

• For the Gaussian distribution an = 2* and = —(3/2) 2* {i — 1). 

• For the uniform distribution an = (2/3)' and ai_(j_i) = -(11/10) (2/3)* {i - 1). 

• For the double-triangular distribution an = 1 and a, = —(17/18) {i — 1). 

Obviously this is no rigorous check from first principles, but if we believe in the regularity and that 
we can at least calculate the first few orders correctly, it adds confidence that no mistake was done at 
higher orders and that the numerical errors are not exceedingly large. A more comprehensive check 
for numerical rounding errors was already mentioned above in this section. 

Table 1: Series for the Bimodal distribution on the d-dimensional hyper-cubic lattice and for x = 



Terms of the series. Xba = 1 + • • • 


+2 


-fx^d' 


+4 x^cP 


+^x^d^ 
+8 x^ d^ 


-fx^d!" 


-40 x^d^ 


+fxU' 
+16 x^d'^ 



2557316 ^5^1 



14175 



+32 X° 



+ ^X^d^ -^x'd"" 
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Table 1: Series for the Bimodal distribution, (continued) 



Terms of the series. 



15891824 ^6^1 



93555 



111488 „6^3 
189 " 



9373372 6 j2 
14175 " 

+64 x^ d^ 



1 190090194848 ™7 jl 
' 42567525 " 

1 645088 ™7 j3 
1 w 

+ 1280 ^7^5 


13641704 ™7 j2 
2079 " 

1 283264 7 M 
' 189 " 

-640 x'^ d^ 


+128 x^d^ 




545049148646 „,8 jl 
127702575 ^ " 

32973784 8^3 

2079 

+ ~ ^8^5 


1 3811431542104 8 ^2 
' 212837625 " 

7834448 ™8 j4 
4725 " 

+ 22016 ^8^6 

+256 d*^ 



2171514982687276 ^9 jl 
8881133625 " 


1 35779921623392 
' 76621545 


X9d2 


32560925165624 9 j3 
127702575 " 


1 3151565216 
' 93555 


x^ d'^ 


31525376 9 j5 
2835 " 


1 1574144 
' 189 


x^ d^ 


1 203392 ^9 ^7 
1 45 X u 


10240 
3 


xU^ 


+512 dp 







202257782879679928 
1856156927625 


x'U' 


3150596158319108 _.10 j2 
7753370625 " 


1 17497504604224 

' 18243225 


d^ 


44042569593584 10 j4 
91216125 " 


1 255487168 
4455 


x^° d^ 


28463552 „10 ^6 
675 " 


1 474880 
' 27 


x^W 


+ 193792 ^10^8 


-7680 


x^U^ 


+1024 



44286591649508625456608 
2143861251106875 

2975087273719088 



,11 



80405325 

14408475958592 
14189175 

270251392 
2025 

I 175104 



d^ 



X 



11 



X 



11 d9 



_ 24516788251206488696 ^11 ^2 



519723939735 

49132009 i:-i94464 
42567525 

7 53069824 
6237 

6347776 
189 

51200 



+ 



;lld6 



X 



11 d« 



X 



11^10 
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Table 1: Series for the Bimodal distribution, (continued) 



Terms of the series. 



+2048 x^^ d 



1 102687986431081211931032 
' 2275791174570375 


^12 Jl 

X a 


98234614240598344870804 
2143861251406875 


^12 j2 

X a 


1631193202018689472 
39978764595 


™12 j3 

X a 


1 132365195242707824 
2170943775 


^12 j4 
X a 


36568729916000 
1702701 


X CI 


1 5544027735104 
' 2837835 


^12 j6 

x u 


1 96053248 
' 297 


^12 M 
X CL 


24126464 
63 


™12 j8 

X a 


1 3447296 
' 63 


^,12 j9 

X u 


1 821248 
1 9 


™12 jlO 

X u 


112640 
3 


12 jll 

X d 


+4096 


12 jl2 

X d 


110468581411293350924457444112 
48076088562799171875 


1 3 jl 


1 172781445528971814087368 
' 28988129795625 


1 3 j2 

x'^^ d^ 


2939279178242203350187328 
510443155096875 


™13 j3 

X a 


1 4778773293081621239776 
' 1856156927625 


™13 M 
X CL 


23565892697470112 
46990125 


™13 j5 

X a 


1 5353769766272 
' 289575 


™13 j6 
X CL 


1 39918090277888 
' 10135125 


^13 j7 

X a, 


1 688011776 
' 693 


™13 j8 
X CL 


179936768 
175 


™13 j9 

X a, 


1 1286144 
+ 21 


™13 jlO 

X ot 


, 1153024 
+ 5 


™13 jll 
X u 


— oiyzu 


™13 j12 

X a 


_i_Qi no 
-|-oiyz 


™13 jl3 
X CI 






305463146085574972582952872664 
33041384503160158125 


14 j1 


1 50334647726118167100558221079016 
' 3028793579456347828125 


14 j2 

x^* 


824015715606745029679283608 
147926426347074375 


^14 j3 

X a 


157096153765758568475862928 
32157918771103125 


^14 j4 
X u 


1 30648160725228437701792 
' 7795859096025 


™14 j5 

X a 


2624211856887202496 
2960377875 


™14 j6 
X CL 


1 43529789374208 
1563705 


T,14 j7 

X a 


1 1037253983072768 
127702575 


™14 j8 
X u 


1 1206767104 
' 385 


^14 j9 
.X CI 


1487067136 
567 


^14 jlO 


5660672 
189 


™14 Jll 
X o 


1 8531968 
1 15 


™14 j12 

X a 


532480 
3 


^14 ^13 

X a 


+iDoo4 


^14 Jl4 

X a 


1 1293826468733333294113597991611697536 
' 3952575621190533915703125 


x^s 


338533911194409672724356315499808 
363455229534761739375 


Xl5d2 


1 3142897456730380314538798434325456 

' 3028793579456347828125 


x^^ 


257903955006958171637687615776 

443779279041223125 




1 5449914292983317242598754368 
32157918771103125 


x^s 


67310752760322121480064 
2998407344625 


x^^d^ 


1 11493150389809720832 
' 23260111875 


x^^d^ 


1 4618407513611776 
' 127702575 


xiSrfS 


1 3526208986206208 
' 212837625 


X^^dP 


1 904493910016 
' 93555 


x^^dio 


90834735104 
14175 


x'^d'' 


93716480 
189 


Xl5dl2 
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Table 1: Series for the Bimodal distribution, (continued) 



Terms of the series. 



_|_ 61898752 ^15 ^13 _ 1146880 ^15 ^14 

+32768 x^^ d}^ 



7 Analysis of the Series 

Our analysis uses the Dlog-Pade method [25] and the methods Ml and M2 [26,27]. Each of these is 
combined with Euler-transformations for improved results. For each series, our main goal is to obtain 
the critical value Xc and the critical exponent 7 which describe the power law divergence, as in 

X^^^A{x^-x)-\l + B{x^-x)^'). (75) 

The scries analysis is done for a fixed dimension at a time. Wc present our results for dimensions 
7 and 8 above the upper critical dimension and for 5 and 4 below it. We also attempted an analysis 
in the physical dimension 3 but the results are not conclusive. 

At the upper critical dimension c/c = 6 the corrections to scaling become logarithmic and there one 
expects the general form 

X^Ax) ~ A{x^ - a;)-^| ln(xc - x)f . (76) 

Instead of Ml and M2, one can apply a modified method to take such corrections into account. This 
was pursued in [8], for the Bimodal distribution, but the authors reported poor convergence already for 
that case. Given that our scries for the other distributions are more problematic, we did not attempt 
a detailed analysis in d = 6. 

It is generally observed in series analysis, that for a given order of expansion, a series behaves better, 
the higher the dimension. That is also the case in the study at hand. Qualitatively it is understood by 
the fast increase of the embedding constants with increasing dimension. Thus a much larger number 
of lattice configurations contributes to the higher dimensional series, allowing it to capture more of the 
underlying Physics. 
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Table 2: Series for the Gaussian distribution on the d-dimensional hyper-cubic lattice and for x = 



Terms of the series. Xea = 1 + • • • 


+ Z X 


a 










a „2 
— D X 


a 


+4 a; 








+21.33333333333333333 


ji 
a 


—24 X 


d 


+8 x^ 


d^ 


—57.9999999999999999 x 


ji 
a 


+93.333333333333333 a; 


j2 

a 


i'7'i A A A A A A A A A A AAA A A 4 

—71.999999999999999 x 


d"' 


+15.999999999999999 a; 


a 










+20.266666666666666 a;° 


ji 
a 


— 119.999999999999999 x° 


a 


1 OKA AAAAAAAAAAAAAAA K 

+359.999999999999999 a;^ 


j3 


— 192 X''' 




+32 x^ 


a° 






+558.4000000000000 x° 


il 

a 


—934.488888888888 x° 


j2 

d^ 


AA-1 AAAAAAAAAAAA R 

—991.999999999999 x° 


d^ 


+1210.666666666666 


74 


^ A A A A A A AAA AAA fS 

—479.999999999999 x° 


d° 


1 O AAAAAAAAAAAA R 

+63.999999999999 x° 


d° 


+1000.02539682539 a; 


a 


; OO^ A O AAAAAAAAAA T 

+2374.39999999999 x 


j2 

d 


+2408.5333333333 X 


d 


-4704.0000000000 a;^ 




+3690.6666666666 a;^ 




-1152.0000000000 a;^ 


d^ 


+128.00000000000 a;^ 












-31435.219047619 a;^ 




+21442.41269841 a;^ 


d2 


-20110.400000000 a;^ 


d^ 


+12771.200000000 x^ 




-18143.999999999x8 


d^ 


+10495.999999999 a;^ 




-2687.999999999 


d' 


+255.9999999999x8 


d8 






-20387.53298060 .x^ 




+294403.1746032 x^ 


d^ 


-95832.1693122 x^ 


d3 


+21440.0000000 x^ 




+61560.888888 x^ 


d^ 


-62080.000000 x^ 


d6 


+28373.333333 x^ 




-6144.0000000 x^ 


d^ 


+512.00000000 x^ 


d^ 



+2051214.7843386 a;i° 


di 


-6192978.8227 xi°d2 


+3181210.20952x1° 


d3 


-1079716.757669x1° 


d* 


-48661.333333 xi°d5 


+257267.199999x1° 


d6 


-195840.00000x1° 


d' 


+73813.33333 xi° d^ 


-13824.00000x1° 


d9 


+1024.000000 xi° 


dio 








-2725463.2041 x" 


di 


+18390206.21 x" d^ 


+2845075.7280 x" 


d3 


+798591.1877 x" 


d4 


-340330.5315 x" d'"^ 


-554547.199x11 


d6 


+961186.1329 x" 


d^ 


-582143.999 xiid8 


+ 186367.999x11 


d^ 


-30719.999x^1 


dio 


+2047.9999 x" rfH 






-52046138.6x12 


di 


+64425096.1 xi^ 


-314627076.1 xi^ 


d^ 


+126856674.2x12 


d* 


-20441160.8 xi2rf5 


-520014.1 xi2 


d6 


-3007795.2 xi^ 


d^ 


+3301421 xi2 d^ 


-1653248x12 


d9 


+459434 xi2 


dio 


-67584.0 xi^d" 


+4096.00 xi2 


dl2 


-2052218007 x" 


dl 


+563481178elxi3 d^ 


-20936942e2 xi^ 


d3 


+846300080 xi^ 


d4 


-15546990elxi3 d^ 


+318945e2xi3 


d6 


+2724596 x" 


d^ 


-1300363elxi3 d^ 


+1063195elxi3 


d^ 


-45281e2 xi^ 


dio 


+1110e3 xi3 (ill 


-147456 xi3 


dl2 


+8192.0x13 


dl3 
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Table 3: Series for the Uniform distribution on the d-dimensional hyper-cubic lattice and for x = 



Terms of the series. Xea = 1 + • • • 


H-U.OOuOOOuOOOuODDODDDO X 


a 










r\ /I o o o o o o o o o o o oo o o oo ™2 


a 


1 ('\ AAAAAAAAAAAAAAAAAAA ™2 

+L). 4444444444444444444 x 


a 






+0.c559 78oc559 78oo597o8 x 


ji 
a 


—0.65185185185185185 a; 


a 


+0.296296296296296296 a; 


d 


+0.08084656084656084 x 


ji 

(X 


+0.373051146384479718 a; 


j2 

a 


—0.65185185185185185 x 




+0.1975c508 64iy75c}08D x 


fit 










—0.89147987814654481 


ji 
a 


+0.84242210464432686 x^ 


72 

a 


+0.49683715461493239 x^ 


j3 

d'' 


—0.57942386831275720 x'' 


a 


+0.13168724279835391 x° 


(T 






+1.0637968983789089 a;° 




— 1.862678750213141 a;° 


l2 

a 


+0.59106251224769744 x° 


d^ 


+0.60287673917303544 a;° 


a 


—0.48285322359396431 a;° 


dr 


+0.08779149519890260 x° 


76 

d° 


+ 1.255336691527167 x 


ji 
a 


— 1.283045789797112 a; 


j2 

a 


—0.625225076837070 X 


d 


+0.3268519171729049 x'^ 




+0.6538376118623032 a;^ 




-0.3862825788751715 a;^ 


d^ 


+0.0585276634659350 x"^ 












-2.764322363362706 x^ 




+5.381812025290449 x^ 


d2 


-1.970220360622241 a;^ 


d^ 


-1.100189850013248 x^ 




+0.06329348749101851 x^ 




+0.6510505802686874 a;^ 


d6 


-0.300442005791800 




+0.03901844231062340 x^ 


rfS 






-9.85782657798808 x^ 




+17.54278833330664 x'^ 


d2 


-9.714423500532811 


d3 


+3.026748425908237x9 




-1.234630286704749 x^ 


d^ 


-0.168708312466884 x^ 


d^ 


+0.608947822994459 x^ 




-0.22890819488898 x^ 


d^ 


+0.0260122948737485 x^ 


d9 



+15.5044273681052 x^^ d^ 
-12.0316553472705 x^^ d^ 
-0.348044505410769 x^^ d^ 
+0.0173415299158332 x^^ d^^ 


-38.7783044835667 x^^ d^ 
+2.84533673882951 x^" d^ 
+0.543557868404698 x^^ d^ 


+33.5641685051471 xi° d3 
-1.14514652844366 xi° d^ 
-0.171681146166746 x^^ d^ 


+ 108.085493498972 x"di 
-48.1504957504608 x" d* 
-0.924877616717460 x^^ d^ 
-0.127171219382772 x^^ d^° 


-232.015426078230 x" d2 
+2.24480229556376 x" d^ 
-0.467778685821369 x" d^ 
+0.0115610199438887 x" d" 


+168.045218530949 xiid3 
+2.83058482344332 x" d^ 
+0.468089181785302 x" d^ 


-142.76066452364x^2 
+231.20946538453x^2^^ 

+2.7869356691639 x^^ d^ 

+0.39186352362185 x^^ rfio 


+418.50354734964 x^^ 
-49.940080322749 x^^ 
-0.64727598233428 x^^ d^ 

-0.093258894214093 x^^ d" 


-458.98712418598 d3 
+0.060628792921663 xi^ d^ 
-0.53174415759476 xi^ d^ 

+0.0077073466292662 xi^ di^ 


-1391.033864136x13 d^ 
+ 1234.36058382x13 d* 
-1.45093243110x13 d^ 
-0.5497466842 x^^ ^jio 

+0.005138231086 x^^ £^13 


+3401.571257928x13 d^ 
-211.0538211516x13 d^ 
+2.63953125141x13 d^ 
+0.32066232142 x^^ d^i 


-3040.190268119x13 d3 
+5.814822817330x13 d^ 
-0.365539203648 x" d^ 
-0.0678246503365 xi3 d^^ 
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Table 4: Series for the Double- Triangular distribution on the d-dimensional hyper-cubic lattice and for 
X = {J/k^Tf. 



Terms of the series. Xea = ! + ••• 


+ i X 


a 










— U. 944444444444444444 x 


^1 
a 


+ 1 a; 


J2 

d 






-1-0.883333333333333333 


ji 
a 


— 1.88888888888888888 a;^ 


a 


+0.999999999999999999 x"' 




-1-0.925282186948853615 x 


ji 
a 


+0.908641975308641977 a; 


j2 

a 


—2.83333333333333333 x 


d^ 


-1-1.00000000000000000 X 


j4 

a 










—6.21639476778365667 x° 


jI 

a 


+7.16815696649029982 x^ 


a 


+1.82592592592592591 x'' 


d^ 


—Aniiiiiiiiiiiiiui x° 
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7.1 Dlog-Pade Analysis 

The Dlog-Pade method is one of the most common methods for the asymptotic analysis of power series. 
One calculates Pade approximants to the logarithmic derivative of the series and obtains estimates for 
the critical value Xc of the expansion variable x (the threshold) and for the critical exponent 7 from 
their real first order poles and the corresponding residues. We also refer to the pole-residue pairs as 
data-points since we often plot them in diagrams of residues versus poles. 

Many series point to singularities other than those representing the physical critical point. They 
are observed in the Dlog-Pade analysis of the original series and, depending on their strength and 
location in the complex plane, hamper convergence of the data points. This effect appears to be 
strongest when an extra singularity is on the negative real axis closer to the origin than the physical 
one. Application of an Euler-transformation into the new variable z = .x„ x/[xn — x)^ with x„ at or 
close to the disturbing singularity, usually improves the behavior of the transformed series. 

For some series, in particular those in higher dimension, we obtain satisfactory results in this 
manner. Data points in the pole-residue plots are high in number and well concentrated along a 
distinct line for each series, examples of which follow below. But for other series, the Dlog-Pade 
method, even in combination with an Euler-transformation, is insufficient for a quantitative analysis. 
So our strategy is in general to use the Dlog-Pade method only to get rough estimates for the critical 
parameters, as a starting point for a detailed analysis with Ml and M2, and to assess the general 
behavior of the series from the number of pole-residue pairs which are obtained. 

7.2 Estimation of Xc and the Critical Exponents Using Ml and M2 

The analysis algorithms Ml and M2 allow the accurate simultaneous determination of the threshold 
Xc, the leading critical exponent 7, and the confluent correction to scaling exponent Ai, assuming the 
asymptotic form 

X{x) - A{x^ - x)-\l + Bix, - (77) 
In Ml, one studies the logarithmic derivative of 

Fix)=jx{x)-{xc-x)^ (78) 

which has a pole at Xq with residue —7 + Ai. For a given trial value of Xc one obtains graphs of Ai 
versus 7 for all Pade approximants of F, and chooses the triplet Xc,7, Ai for which best convergence 
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of the different approximants results [27]. 

In the M2 method one first transforms the series in x into series in the variable y = l — {l — x/xc) 
and then takes Pade approximants to 



which should converge to —7. Here one plots graphs of 7 versus the input Ai for different trial values 
of Xc and again chooses the triplet Xc,7, Ai with the best convergence of all Pade approximants. For 
both methods it is advisable to perform first the usual Dlog-Pade analysis, to get rough estimates of 
Xc and 7 which one uses as starting points for the detailed analysis with Ml and M2. The effectiveness 
and preciseness of these series analysis methods has been demonstrated in several papers [26,28-30]. 

In Ml we vary the trial-Xc until the curves from the high order Pade approximants enter fairly 
symmetrically from both sides and the best convergence is obtained. This Xc and the corresponding 
7 are taken as the temporary best estimates for that series, with temporary error estimates from the 
nearest trial-Xc's, whose plots show poorer convergence. In many cases Ml proves to be quite sensitive 
to small changes in the trial-Xc, and the degree of convergence usually looks very convincing. Away 
from the best Xc, convergence degrades quickly, the picture becomes non-symmetric and at the same 
time the area of convergence shifts to lower or higher values of 7. We show examples of such plots 
in Sec. 17.41 In M2 we vary Xc and look for best convergence of the Pade approximant curves while 
they cross each other with a small negative slope. Compared to Ml, the M2-plots are often much 
less decisive. A good convergence region sustains over a wider range, where again the change in Xc is 
accompanied by a shift in the corresponding 7. 

In the end we determine an overall estimate for Xc, which is consistent with the estimates from 
both Ml and M2. These numerical results are presented in the tables of Sec. 17.41 In the tables we 
also include rough estimates for Ai. We comment that the Euler transformation is known to produce 
analytic correction terms even if not present originally. When the leading correction exponent is larger 
than 1, as seems to be the case for some of our series, these 'artificial' corrections will show up in Ml 
and M2 [31], and hence our Ai estimates are mainly included for reference and should not be trusted 
as the real physical values. 




(79) 
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7.3 Sensitivity to the Parameter of the Euler Transformation 

Our analysis relies in a large part on the use of Euler transformations to increase the number of useful 
Pade approximants and to improve their convergence. The technique is well established and has been 
used with success [32], but nevertheless we find it worthwhile to check, to what degree our results are 
sensitive to the precise choice of the parameter Xn, the value of x that is mapped to infinity by the 
transformation. We first choose very close to the negative singularity, as indicated by the Dlog-Pade 
analysis of the original series. We then vary this x-a over a considerable range of typically 20%, and 
compare the results. We observe that a variation of x^ does move the data points or curves obtained 
from individual Fade approximants, but that the average (in Dlog-Pade plots) and the convergence 
region (in Ml plots) stay fixed to a very good accuracy, when compared to the error bounds given 
by the analysis technique itself. We thus exclude that our results are artifacts of the applied Euler 
transformations. 

7.4 Exphcit Results from the Analysis 

The numerical results for dimension 8 are summarized in Tab. [51 In this dimension, even without an 
Euler transformation, the Dlog-Pade analysis gives convincing results for all the distributions: Bimodal, 
double-triangular, uniform and Gaussian. 

During the analysis we prepared a large number of plots of which we can only present a few to 
illustrate the process. The distribution of the pole-estimates indicates a negative real pole for the 
asymptotic function, somewhat weaker than the positive one (e.g. Fig. Both poles are at the same 
distance from the origin. The convergence of pole-residue pairs improves upon Euler transformation. 

The numerical results for the critical exponents are calculated as averages over estimates from 
high-order Pade approximants. These include data from the untransformed series and from the series 
transformed with three different values of Xn- In the tables these are the entries without an estimate 
for Ai. All values are slightly larger than the mean field value of 7 = 1 (Fig. [l]). This deviation is 
understood on theoretical grounds as being due to corrections to the leading singular behavior. Also 
the data indicate that the exponent estimate may further approach 1 for longer series, since generally 
residues decrease in value as the approximant-order increases while remaining greater than one (Fig. [2)) . 
When taking into account the correction exponent with M1/M2 the deviation also decreases. The 
results from Ml and M2 are shown in the table separately for different values of Xn, and generally 
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show very good agreement. We observe the possibihty of a systematic shift to 7 = 1.060(12) in case 
one chooses a different region of best convergence. The values for Ai are included for reference only, 
due to the reason mentioned before. Although the absolute value of the exponent-estimate is larger 
than 1, we find numerical agreement of the results for all the tested distributions. 

Table 5: Results for dimension d = 8 from the analysis with Dlog-Pade, Ml and M2. The first line for 
each distribution shows the result from the Dlog-Pade analysis in which Euler transformations with 
different values were used. The remaining lines show the results from Ml in combination with M2, 
separately for several values of Xn- 



Distribution 


Parameter Threshold Exponent Correction- 
Xn a^c 7 Exponent 


Bimodal 


several 0.072 1.05(1) n/a 
-0.084 0.07331(3) 1.046(9) 1.4-1.7 
-0.073 0.07331(3) 1.046(9) 1.4-1.7 
-0.056 0.07332(3) 1.047(12) 1.3-1.5 


Gaussian 


several 0.080 1.068(20) n/a 
-0.084 0.08030(3) 1.048(9) 1.3-1.5 
-0.070 0.08029(3) 1.047(12) 1.3-1.5 
-0.056 0.08030(3) 1.048(9) 1.3-1.5 


Double- 
Triangular 


several 0.148 1.072(22) n/a 
-0.168 0.14895(3) 1.048(9) 1.3-1.7 
-0.140 0.14895(3) 1.048(9) 1.3-1.5 
-0.112 0.14898(9) 1.048(9) 1.3-1.5 


Uniform 


several 0.228 1.069(25) n/a 
-0.252 0.22852(9) 1.048(6) 1.3-1.5 
-0.210 0.22848(6) 1.048(6) 1.3-1.5 
-0.168 0.22854(9) 1.048(6) 1.3-1.5 



The qualitative behavior in dimension 7 agrees with that for d = 8, although the exponent estimates 
are slightly farther away from 1 (see Tab. IHl). The critical threshold Xc for each series is larger than 
in dimension 8. Again we observe a negative pole of comparable strength and distance from the 
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Figure 3: Position of the Poles from all Pade approximants to the logarithmic derivative of the un- 
transformed series of Xea ^o^' uniform distribution in c? = 8. The upper half shows the location 
of the poles in the complex plane. The lower part is a histogram of the poles along the real axis, in 
which the two peaks are of relevance. The one on the positive axis represents the physical singularity 
we want to characterize, the second one can interfere with analysis and may be mapped away using 
an Euler transformation. 
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Figure 4: Pole-residue plot from a Dlog Pade analysis of the Xea"^^^^^^ uniform distribution in 

d = 8. We use Euler transformations with several values p„. The main part gives an overview over a 
large region including almost all data points. The histograms on the axes show how these points are 
distributed. The inset is an enlarged view of the small region with the highest concentration of points 
as indicated by the box. For comparison we calculate the average and standard deviation c7„ from the 
points in the boxed area. 
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Figure 5: Estimates for the leading critical exponent 7 from the Xea"^^'^^^^ double-triangular 
distribution in d = 8. Series transformed with different values for Xn, indicated by different symbols, 
contribute to the plot. Each point corresponds to a particular transformation and Pade approximant 
[L/M], where L and M relate to the number on the x-axis as given in the label. The different symbols 
refer to different parameters for the Euler transformation: Without (+), Xn = —0.140 (x), Xn = -0.168 
(o), and Xn = —0.112 (*). Estimates are calculated as averages over points in the convergence region. 
The points that are included in the average are those within the boxed area. The resulting numbers, 
in this case of 7 = 1.07 it 0.02, enter table |S1 
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Figure 6: Pole-residue plot from a Dlog Pade analysis of the untransformed Xea"^^^^^^ Uniform 
distribution in d = 7. See the caption of Fig. 0] for further explanation. 
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origin as the physical one, but nevertheless, even without an Euler transformation the series give 
consistent results. To illustrate this we show in Fig. |21 a pole-residue plot of the untransformed series 
for the uniform distribution and also include in Tab. IHl the estimates from the Dlog-Pade analysis 
without transformation. In our Ml/M2-analysis we again observe the possibility of a systematic shift 
to 7 = 1.120(15) in case one chooses a different region of best convergence. 

In dimension 5 the negative pole is closer to the origin than the positive one, and convergence 
degrades. Still the data points line up properly. To improve their convergence and to get higher 
numerical accuracy we apply Euler transformations with several (e.g. Fig.[7| and the final estimates 
(Tab.[7I) are obtained the same way as described before. Looking at Fig.|HJ the main line of data-points 
for the Gaussian distribution still increases with the order of the Pade approximants. Therefore, the 
value measured by this method is probably smaller than what a longer series would show. In summary, 
all studied distributions agree at d = 5 on a common exponent within their error margins. 

In Figures |H1 to we show, for the case of the Gaussian distribution, plots as they are typically 
obtained from Ml. Each curve in a plot comes from a different Pade approximant as the legend shows. 
In all figures we clearly see a region where the lines converge, and since all figures show the same range 
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Table 6: Results for dimension d = 7 from the analysis with Dlog-Pade, Ml and M2. The first line for 
each distribution shows the result from the Dlog-Pade analysis in which no Euler transformations was 
applied. The remaining lines show the results from Ml in combination with M2, separately for several 
values of x^. 



Distribution 


Parameter Threshold Exponent Correction- 
Xn Xc 7 Exponent 


Bimodal 


none 0.088 1.14(3) n/a 
-0.082 0.08731(9) 1.105(15) 1.5-1.7 
-0.078 0.08732(9) 1.105(12) 1.3-1.7 
-0.051 0.08738(6) 1.110(15) 1.3-1.4 


Gaussian 


none 0.097(1) 1.14(3) n/a 
-0.071 0.09710(6) 1.107(9) 1.3-1.4 
-0.062 0.09712(9) 1.108(12) ^1.3 


Double- 
Triangular 


none 0.1784(4) 1.14(3) n/a 
-0.265 0.17790(9) 1.108(9) 1.3-1.6 
-0.156 0.17790(9) 1.108(9) 1.3-1.5 
-0.112 0.17799(9) 1.110(6) ^1.3 


Uniform 


none 0.2745(5) 1.13(3) n/a 
-0.315 0.27402(9) 1.108(9) 1.3-1.5 
-0.234 0.27399(9) 1.108(6) 1.3-1.5 
-0.225 0.27399(9) 1.108(6) 1.3-1.5 
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Figure 7: Pole-residue plot from a Dlog Pade analysis of the Xea"^^'^^^^ Gaussian distribution 

in d = 5. We use Euler transformations with several values p„. See the caption of Fig. 0]for further 
explanation. 
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in the 7-Ai-plane it is also easy to see that the convergence region shifts around. Figures Ufa) and 
(b) differ in the trial x as input parameter. While for x = 0.177 convergence is quite good, is becomes 
better for x = 0.179 (in fact best among our trial values). At this value of x the shape of the curves is 
also symmetric and they switch over to one side for larger x (the opposite side, when compared to the 
smaller x) which we find to be a characteristic feature of the point of best convergence. Figure [TUl shows 
the corresponding plot for the untransformed series. Here x = 0.179 is also near the characteristic point 
of symmetric curves. Still convergence is not as good as in Fig. IHtb) and the center is shifted to a 
slightly larger value of 7. The example, although not the most common case, also illustrates that one 
must not rely on analysis with either Ml or M2 alone. While Ml gives an estimate of 7 = 2 or higher, 
M2 (Fig. Ill|) points to a lower value of roughly 1.82 and Ai above 1, or 7 ~ 1.95 with Ai below 1 
and poorer convergence. Our estimates in the tables always result from using Ml together with M2. 
Fig. ^1 is another example, showing a plot from Ml-analysis for the Uniform distribution near the 
symmetry point of best convergence. 

As we decrease the dimension further to d = 4, the analysis becomes increasingly difficult. One 
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Table 7: Results for dimension d = 5 from the analysis with Dlog-Pade, Ml and M2. See caption of 



Distribution 


Parameter Threshold Exponent Correction- 
Xn Xc 7 Exponent 


Bimodal 


several 0.154 1.91(10) n/a 
-0.120 0.154(3) 1.95(15) 1.1-1.3 
-0.100 0.154(3) 1.95(15) «1.0 


Gaussian 


several 0.174 (1-67(8)) n/a 
-0.096 0.176(3) 1.70(15) 0.8-1.0 
-0.080 0.177(3) 1.75(15) 0.8-1.0 
-0.064 0.177(3) 1.75(15) 0.8-1.0 


Double- 
Triangular 


several 0.312 1.81(7) n/a 
-0.240 0.312(6) 1.80(15) 0.9-1.0 
-0.200 0.312(6) 1.80(15) 0.9-1.0 


Uniform 


several 0.484 1.72(6) n/a 
-0.348 0.484(6) 1.70(15) 1.0-1.2 
-0.290 0.487(6) 1.70(15) 0.8-1.0 
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Figure 8: Estimates for the leading critical exponent 7 from the XsA'S^'^i^^ Gaussian distribution 

in d = 5. See the caption of Fig. [3 for details. The distribution of the data-points indicates that the 
7-estimate may still increase with longer series. The different symbols refer to different parameters 
for the Euler transformation: Without (+), Xn = —0.080 (x), x„ = —0.096 (o), and x„ = —0.064 (*). 
Averaging over the points inside the small box leads to an estimante of 7 = 1.68 it 0.08, which enters 
table m 
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Figure 9: Ml analysis of the x^^-series for the Gaussian distribution in d = 5. An Euler transformation 
with Xn = 0.08 was apphed. Plot (a) was obtained for a trial x of 0.177, where convergence is visible 
but lies below the optimal convergence point. In (b) the trial x = 0.179 is near or at the point of 
optimal convergence. 

(a) 



2.2 r 



1.9 - 



a 
o 

X 



1.7- 



1.6 - 



^ [4/6] 
O [4/7] 

^- [5/5] 
* [5/6] 

[6/4] 
+ [6/5] 
[7/4] 




1.5 - 



1.4 L 



0.6 0.7 0.8 0.£ 



1 1.1 1.2 1.3 

Correction Exponent Ai 



1.5 



1.6 



(b) 




^ [4/6] 

O [4/7] 

-^- [5/5] 

* [5/6] 

-e- [6/4] 

+ [6/5] 
[7/4] 



1 1.1 1.2 1.3 

Correction Exponent Ai 



1.5 



1.6 



39 



Figure 10: Ml analysis of the untransformed Xea"^^^^®^ Gaussian distribution in d = 5. The 

trial X = 0.179 is near or at the point of optimal convergence. Compared to Fig|U[b), with Euler- 
transformation, convergence is less sharp and shifted to a slightly larger 7. Strong fluctuations of the 

Pade-lines are visible in the lower part of the plot. 
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Figure 12: Ml analysis of the XsA'S^'^i^s Uniform distribution in d = 5, after an Euler- 

transformation with = 0.29. x = 0.488 is near or at the the point of optimal convergence. 
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reason is the negative pole on the real axis, which is very strong and apparent for all the series. 
Without an Euler transformation, the Dlog-Pade analysis does not show anything conclusive. For the 
transformed series, a larger number of data points lines up well in the pole-residue plots, but they 
are still not well converged. The series for the Gaussian distribution is somewhat exceptional here: 
Exponent estimates converge well with increasing approximant-order, but they approach a value of 

7 = 3.1 ±0.1 which is much higher than what we obtain in the other cases and with the other methods. 
If this is not simply an artifact of the transformation, it must be attributable to the correction-term 
to scaling, which becomes increasingly important at lower dimension. 

With M1/M2 applied to the transformed series we are able to obtain estimates for Xc and 7, 
although the error margins are quite large. Indeed, the value we obtain for the correction exponent 
Ai is much larger than 1 (and larger than in the higher dimensions). Again we find that our numbers 
agree with a common exponent 7 for all the distribution functions. 

8 Conclusions 

Figure summarizes our numerical estimates for the leading critical exponent 7 of the Edwards 
Anderson susceptibility Xea different dimensions. For each dimension we show the 4 values 

obtained for the different random distributions. The mean values and error bars are as shown in 
Tables El to 8. 

The error bars for dimensions 7 and 8 are too small to be visible in the plot. We observe once 
more that in these dimensions the estimates are close to, but still larger than the expected mean-field 
value of 1. As we argued, the deviation is likely caused by correction to scaling terms and the fact that 
we work with relatively short series. Our Dlog-Pade analysis suggests smaller estimates with longer 
series, and when accounting for the first correction term using the methods Ml and M2 the deviation 
indeed decreases, but higher order corrections cannot be excluded. The larger deviation for d = 7 is 
consistent with the observation that the correction terms are more dominant in lower dimensions. 

As expected, the estimates we obtain for dimensions 4 and 5 are pronouncedly different from the 
mean field value. Within each dimension the estimates agree on a common value for all the random 
distributions we study, which is roughly 7 = 2.4 it 0.2 in d = 4 and 7 = 1.8 it 0.2 in d = 5. Thus our 
data do not indicate that the random distribution for the quenched-in disorder splits the spin glass 
model into many universality classes nor that the model behaves in that respect differently than other 
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Table 8: Results for dimension d = A from the analysis with Dlog-Pade, Ml and M2. The table 
lines without estimates for Ai stem from the Dlog-Pade analysis in which Euler transformations with 
different values Xn were used. The remaining lines show the results from Ml in combination with M2, 
separately for several values of x^. 



Distribution 


Parameter Threshold Exponent Correction- 
Xn Xc 7 Exponent 


Bimodal 


-0.144 0.26(2) 2.5(3) 1.5-1.6 
-0.120 0.26(2) 2.5(3) 1.5-1.6 


Gaussian 


several 0.31(2) 3.1(1) n/a 
-0.108 0.312(4) 2.3(1) 1.3-1.4 
-0.090 0.314(4) 2.3(1) 1.3-1.4 
-0.072 0.314(4) 2.3(1) 1.3-1.4 


Double- 
Triangular 


several 0.52(8) 2.8(8) n/a 
-0.276 0.54(2) 2.5(2) «il.5 
-0.230 0.54(2) 2.5(2) «1.5 


Uniform 


-0.396 0.84(2) 2.5(1) 1.3-1.4 
-0.330 0.83(2) 2.4(1) 1.3-1.4 
-0.264 0.84(2) 2.4(1) 1.3-1.4 
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common thermodynamic models, in contrast to claims from Refs. [14-20]. Instead we find confirmation 
for the established picture, that the space dimension creates universality classes and that the leading 
critical exponent is a universal quantity [33]. 

Most of the simulations of the Ising spin glass have been done in dimension 3, in which our series 
do not perform well. For a direct comparison we are thus limited to the sparse results for c/ = 4 
from [16, 18]. Our estimates for the critical temperature Tc agree rather well with those by Bernardi 
and Campbell. The compared values are Tc = 1.96 (vs. 1.99 it 0.01) for the bimodal distribution, 
Tc = 1.88 (vs. l.DlibO.Ol) for the uniform distribution, and Tc = 1.79 (vs. 1.77ib0.01) for the Gaussian 
distribution, where the uncertainty in our estimates is also roughly 1 in the last digit, from fluctuations 
and from possible additional systematic shifts, due to scaling corrections. We confirm a slight decrease 
of Tc with increasing kurtosis of the random distributions, which is defined as the ratio of moments 
R = M4/M2. The kurtosis values are: Bimodal 1, double-triangular 4/3, uniform 9/5, Gaussian 3 and 
for the exponential distribution 6. Bernardi and Campbell have calculated the exponent r/, while we 
have values for 7, so we currently lack a third exponent, such as v, for a direct comparison. However, 
the discrepancy in the general universal vs. non-universal behavior remains. 

Some authors [34-36] have stressed the importance of corrections to finite size scaling (FSS). In 
taking these into account, they do not find violated universality. We cannot assess the quality of 
the simulations that were done by Bernardi, Campbell [14-16, 19] and others, but, from the data in 
the papers and later citations, we are led to speculate that neglected corrections to FSS have caused 
systematic errors in the exponent estimates from simulations. We would further like to stress the 
general statement made by other authors, that the characteristic features of the spin glass with its 
quenched-in disorder creates enormous problems for simulations. Series expansion techniques appear 
here particularly suitable since the configurational average over the randomness is handled exactly 
within their framework and own limitations. 
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Figure 13: Estimates for the critical exponent 7 grouped by space dimension. The different values for 
each dimension are obtained from the 4 probability distributions. The estimates appear to agree on 
a common exponent 7 within but separately in each dimension. Above the critical dimension = 6 
the values are close to 1. 
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